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Abstract 


Mesons are made of quark-antiquark pairs held together by the strong force. 
The one channel spectator, Dirac, and Salpeter equations can each be used 
to model this pairing. We look at cases where the relativistic kernel of these 
equations corresponds to a time-like vector exchange, a scalar exchange, or a 
linear combination of the two. Since the model used in this paper describes 
mesons which cannot decay physically, the equations must describe stable 
states. We find that this requirement is not always satisfied, and give a 
complete discussion of the conditions under which the various equations give 
unphysical, unstable solutions. 
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FIG. 1. Example of electroweak decay of the meson. 

I. INTRODUCTION 
A. Background 

In the simplest models, mesons are bound states of a valance quark-antiquark pair con¬ 
fined by the strong force. Even for such a simple case a covariant model is needed when the 
mesons are composed of light quarks with high momentum components. However, covariant 
models require knowledge of the Lorentz structure of the conhning interaction, and it turns 
out that some choices of Lorentz structure for some equations will produce mesons which 
decay. When no mechanism for decay has been included in the model (which will be the sit¬ 
uation for the cases discussed in this paper), this is a sign that the solutions are unphysical. 
In may be acceptable for an equation to produce unstable (i.e. unphysical) solutions /it if 
these solutions are conhned to a region of the spectrum which can be precisely characterized 
and systematically ignored, but if this is not possible equations which produce such unphys¬ 
ical solutions are unsatisfactory. In this paper we study conhning potentials with scalar and 
time-like vector exchanges, and hnd that the stability of such interactions depends on the 
kind of relativistic equation used for the description of the interaction. 

This is not the hrst time that the stability of covariant models of conhnement has been 
addressed. Several papers have been written on this topic, some with contradictory con¬ 
clusions. Two examples which illustrate this are papers titled An exact argument against 
an effective vector exchange for the confining quark-antiquark potential [|^, and Evidence 
against a scalar confining potential in QCD p|. If both papers are correct, this would in¬ 
dicate that, at best, the Lorentz structure for the potential is more complex than a simple 
scalar or vector exchange. 

Our research into the question of stability was motivated by the paper of Parramore and 
Piekarewicz , which found that the Salpeter equation was stable when the vector strength 
exceeded the scalar strength. This seemed counter intuitive to us, since it is well known that, 
because of the famous Klein paradox, the Dirac equation is stable only when the potential 
is predominately scalar. Their result also contradicted the work of another group [Q who 
found that the Salpeter equation was stable for a pure scalar conhning interaction, provided 
the quark mass was sufficiently large. 
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FIG. 2. Example of strong decay of the p"*" meson. 

B. Physical and unphysical instabilities 

We begin the discussion by making a distinction between instabilities which are physical 
and those which are unphysical. Real mesons have a finite lifetime and can decay either 
through the strong interaction or the electroweak interaction. For example, the can 
decay into a photon and a through the electroweak interaction shown in Fig. |l|. It 
can also decay into a tt’*' and 7r° via the strong interaction, as shown in Fig. In this 
paper we describe mesons which are isolated from external influences (including vacuum 
fluxuations), and use an equation which excludes the electroweak interaction and does not 
include any mechanism for the production of quark-antiquark pairs. Hence both of these 
decay mechanisms are excluded from the theory and thus the mesons described by our 
equations cannot decay physically. Therefore any instability emerging from these equations 
will be unphysical, and a sign that the equations are describing unacceptable states. 

C. Unphysical instabilities — an example 

The Dirac equation for a linear combination of a scalar and vector conhning potential 
provides a familiar example of the kind of unphysical instabilities we are discussing. Consider 
the Dirac equation for the linear conhning potential V{r) = ar {(1 — y) + 2/7°} 

EB'y^(p{r) = {m + V{r) P-'y-W) (p{r) (1.1) 

where a and the vector strength y are both constants. The solutions of this equation have 
both positive and negative binding energy eigenvalues Eb- If the system described by this 
equation could interact with the outside world (e.g. absorb or emit photons), the positive 
energy states could decay to negative energy states (unless all of the negative energy states 
were occupied as in hole theory). However, we have assumed that there is no coupling to 
the outside world, and hence this equation should describe a stable system, even if some of 
the binding energy eigenvalues are inhnitely large and negative. However, it is well known 
that the Dirac equation does not give stable solutions for all values of the vector strength y 
and we review this result now. 

The nature of the solutions to the Dirac equation can be studied by looking at the 
expectation value of U = m + V. The form of this expectation value, which describes how 
the wave function behaves, is 
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FIG. 3. Sketch of the solution to the Dirac equation for the scalar case, where a > 0 and y = 0- 



FIG. 4. Sketch of the solution to the Dirac equation for the vector case, where a > 0 and y = 1. 



m + ar 

—m — ar{l — 2y) 


positive energy 
negative energy 


( 1 . 2 ) 


where the positive energy expectation valne is a matrix element involving w-type positive 
energy spinors, (t/)+ = uUu, and the negative energy expectation value is a matrix element 
involving n-type negative energy spinors, (t/)- = vUv. The result (|1.2|) comes from the 
matrix elements 


uu = 1= —vv 

= 1= v'j^v , (1.3) 

which hold when the total momentum p = 0. When Eq. ( 0 ) is sketched for pure scalar 
(y = 0) or pure vector (y = 1) cases, Fig. || and Fig. | are produced, respectively. The 
resulting wave functions for a particle with energy E are sketched on the hgures, along with 
the form of (U) which produces it. 

To understand these results, hrst neglect the coupling between positive and negative 
energy states. Then the positive energy states move under the influence of the potential {U) + 
and the negative energy states under the influence of {U)-. For the scalar case (y = 0), the 
choice cr > 0 produces conhnement for both positive and negative energy states. Coupling 
the two solutions does not change this picture signihcantly, and the exact solution is a total 
wave function which drops to zero at large distances. This means that both positive and 
negative energy solutions describe particles permanently conhned around the point r = 0. 
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Next look at the vector case {y = 1), and begin again by neglecting the conpling between 
the positive and negative energy states. In this case, however, either the positive or negative 
energy state is always nnconhned. For the example shown in Fig. cr > 0 and the positive 
energy states are conhned and the negative energy states are not. Inclnding conpling between 
the positive and negative energy states mixes the two states, and the wave fnnction for the 
exact positive energy solntion acqnires a component with a “tail” which oscillates to inhnity, 
signaling deconhnement. The effect of the conpling is to prodnce an effective potential 
composed of two regions separated by a hnite potential barrier throngh which the qnark can 
tnnnel. Once it is free of the potential barrier it can propagate endlessly throngh space, thns 
becoming a free qnark. In this case, the exact conpled solntions do not conhne either the 
positive or negative energy states, and the bonnd state is nnstable. This example, known 
as the Klein paradox |p, is one of the nnphysical instabilities we are trying to avoid. 


D. Requirements for stability 

A relativisitic eqnation with a conhning kernel with a given Lorentz strnctnre will have 
stable, physical solntions only if the following fonr conditions are satished: 

(1) the binding energy mnst be real; 

(2) the energy eigenvalnes mnst be independent of the nnmerical approximations nsed 
to obtain them; 

(3) nnphysical solntions, if there are any, mnst be conhned to an identihable part of the 
spectrnm clearly separated from the physical solntions; and 

(4) the solntions mnst have the correct strnctnre in coordinate or momentnm space. 

We will discnss each of these conditions in tnrn. 

Condition 1 - real energies. Any eigenstate wave fnnction which describes a meson in 
momentnm space, can be written 

= 4‘(p)e-“^‘, (1.4) 

where E = The discnssion is simplihed if the particle is chosen to be at rest, 

P = 0. Then, if p is complex, y = yoEiV /2, the absolnte sqnare of the meson wave fnnction 
is 


l^(p,^)r = (1-5) 

As time increases, this goes exponentially either to zero or to inhnity, showing that the state 
is nnstable. 

Condition 2 - numerical stability. The diherent relativistic eqnations will be solved 
nnmerically in Sec. IV nsing spline fnnctions to model the wavefnnctions in momentnm 
space. (A description of the properties of the spline fnnctions is given in Appendix A.) 
So long as enongh spline fnnctions are used to model the system, the energy of the lower 
lying stable states will not vary much as the spline rank is increased. However, if the state is 
unstable it is part of a continuous spectrum and the energies obtained from the “eigenvalue” 
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equation only represent a discrete approximation to this continnous spectrum. They will 
vary strongly with the number of splines, mnch as the location of the nth point in the 
interval [0,1] will vary strongly with the number of intervals N into which the the line 
segment is divided. This dependence of an energy level on spline rank is one of the most 
obvious symptoms of instability. 

Condition 3 - isolation of instabilities. In some cases we hnd that, following the second 
criteria, the positive energy states are stable and the negative energy states are unstable. 
This may be acceptable for a phenomenology, where the negative energy states can be 
rejected as nnphysical from the start. However, in some cases these unstable negative energy 
states become positive as the spline number increases, and they can become so positive that 
they cross the gap separating the negative and positive energy states, enter the positive 
energy spectrnm, and mix with states which would otherwise be stable. In this case the 
distinction between (stable) positive energy states and (unstable) negative energy states 
becomes blurred, and we cannot rely on the predictions of the equation. 

Condition 4 - correct structure. Even if the mass is real, the state might not be confined 
in a hnite region of coordinate space (as in the Dirac example ontlined above). If the state 
is conhned, its coordinate space wave fnnction will approach zero as r —cx) faster than an 
exponential. It can be shown that the momentnm space wave fnnction resulting from snch 
a state will also fall off at p —> cx) faster than an exponential, and that the nnmber of nodes 
will correspond to the level of the state. It is easy to distinguish such behavior from that 
of an unconhned state, which is neither localized in coordinate nor momentum space, and 
which has many nodes not related to the level of the state. We can use the Dirac wave 
functions for comparison, since we know that they are stable for scalar conhnement and 
unstable for vector confinement. Examples of both types of states will be given in Sec. IV. 

In the following sections, these stability conditions will sometimes be referred to by 
nnmber, as we will see that a snccessful phenomenology reqnires that all of them be satished. 


E. Summary and Outline 

In snmmary, the stability conditions are: (1) the eigenvalnes of the system must be real; 
(2) the eigenvalues cannot vary with the spline rank; (3) the positive energy states must 
always be greater than any unstable negative energy states; and (4) the wave fnnctions must 
have the appropriate structure for that specihc state. 

In Sec. II specihc forms are given for the Dirac, Salpeter, and one channel spectator 
(denoted ICS) relativistic equations. Then in Sec. Ill these three equations are stndied 
using an approximation technique which gives insight into the origin of the instabilities, and 
the estimated masses of stable states are compared to the exact numerical solutions presented 
in Sec. IV. The three equations are solved numerically in Sec. IV using spline functions for 
a quasirelativistic conhning potential. The actual equations used in the computer code and 
the properties of spline fnnctions are given in Appendix A. Finally, conclusions are given in 
Sec. V. 
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FIG. 5. Feynman diagram for the meson bound state vertex function. The kernel, or potential, 
is denoted by V. 


II. THE RELATIVISTIC EQUATIONS 

In this section we define the one-channel spectator (ICS) equation obtained by conhning 
the heavier particle 1 (assumed to be the quark) to its positive energy mass shell, hxing 
the ko integration. Then we show that these equations reduce to the Dirac equation for the 
lighter particle (particle 2) in the limit when the mass of the heaver particle mi —> oo. We 
conclude by hnding a helicity representation for the ICS and for the Salpeter equation. 


A. Dirac form for the one-channel spectator equation 


The Feynman diagram for the bound state meson vertex is shown in Fig. Particle 1 is 
the quark, particle 2 the antiquark, and 0 is a matrix in Dirac space which describes how 
the conhning force couples to the quark or antiquark. It can be a scalar, 1, or the time 
component of a four-vector, 7 °. The kernel V contains the momentum dependent structure 
of the conhning potential. The equations are derived in the center of mass rest frame, 
P = (/i, 0 ). Later, the quark will be placed on shell, thus producing the single channel 
equation. The four momenta used in the diagram are 

Pi= P + ^P P2= P-\P 

ki = k + ip ^ 2 = k — ip. ( 2 . 1 ) 

The vector k is the average internal momentum and vector p is the average external mo¬ 
mentum of the quark-antiquark pair: 


P=\{.Pl+P2) P = PI-P2- (2.2) 

With this notation, the Bethe-Salpeter equation for the bound state vertex function for 
the meson is 



V{p, k) 0 


mi+ Jki 
ml — ki 


m 


m2+ M ^ 
ml-kl 


(2.3) 


The two fermion propagators have poles in the complex ko plane; these four poles are 
shown Fig. |^. Factoring the denominators of the propagators 
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FIG. 6. This figure shows the position of the four poles associated with the four propagators 
Gf in the bound state equations. 


the poles are at 


pole 1 


pole 3 

(Gr)-' 

pole 2 


pole 4 

(G 2 -)-^ 



= Ek^ - {ko + i/i) - ze = 0 

= Eki + {ko + |/r) — ze = 0 

= Ek^ - {ko - |/z) - ze = 0 

= Ek^ + {ko - |/i) - ze = 0 


ko = Ek^ - i/i - ze 
ko = -Eki -\ii + ie 
ko = Ek^ + ln-ie 
ko = ~Ek2 + \k' + 


(2.4) 


(2.5) 


To place particle 1 on the positive energy mass shell, the ko integration is closed in the lower 
half plane and only the residue from pole 1 is kept. This gives the following equation: 


r(p) 


d?k 

‘2Ek,{27i)^ 


V{p,k) 0 (mi+ M) T{k) 


m2+ h ^ 
ml -kl~' 


( 2 . 6 ) 


where now ki = {Ek^, k) and ^2 = {Ek^ —/z, k). This is one form of the one channel spectator 
equation. 

Next, we recall that the projection operator can be written [Q as a sum over on-shell u 
spinors (to be dehned below) 


^1+ jki = M(k, A')M(k, A'). 

A' 


(2.7) 


Therefore, if we dehne the relativistic meson wave function by 

^{k. A) = -^= zJ(k, A) T{k) 

flir) - /pc 


2K 


ki 


( 2 . 8 ) 


2 '^2 


then Eq. 


becomes 


T(p, A) (m2- h) = - J 


d^k V{p, k) 


(2ir)=> y 


Y,eti:(.p,k)^{k,X)e, 


(2.9) 


where 














(p, *:) = ti(p, A) e !i(k, A'). 


( 2 . 10 ) 


Equation (^) is the Dirac form of the one-channel spectator equation. Later we will 
reduce this equation further, but as written in p.9|) it looks very much like a Dirac equation 
for the light antiparticle (particle 2) moving under the influence of a effective potential 
which depends on the spin of the heavy quark. To make this comparison more familiar, we 
convert the equation into the usual form by taking the transpose and multiplying by the 
Dirac charge conjugation matrix C (see Ref. 0). This gives 


{m2+ h) ^(P, A) = -0 ^ 

A' ' 


dr’k 

(^ 


T(fc,A') 


V{p.k) 

4:Epi 


0 AA+(P,fc), 


( 2 . 11 ) 


where 


'i'(p,A) = C>typ,A) 


e = ce‘c 


T/O—1 


( 2 . 12 ) 


With the exceptions of the spin dependence of the source, expressed through the factor 
(p, /c), and the fact that the effective potential does not depend solely on q = p — k, 
the difference of the three momenta, Eq. ( |2.11|) looks like the familiar Dirac equation for a 
particle with four momentum equal to —p 2 (as expected from the charge conjugate state). 

We will now calculate the matrix element 0’''+ and show that Eq. ( |2.11| ) does indeed 
reduce to a Dirac equation in the limit mi —> cx). In spin space the u spinor, as defined in 
Ref. [0, is 


/ 


m(p, s) = (Ep + m )2 


11 


a-p 


\ 




(2.13) 


\ Ep + m J 


which contains the operator cr ■ p. It is convenient to work in helicity space, where cr-px^ = 
2A|p|x^. The u spinor in helicity space is therefore 


M+(p, Xj) = m(p, Xj) = Nj 


Pj 


( 1 


\ 




(2.14) 


where we have introduced the notation with p-spin = -|-1 for the positive energy solutions 
(u) and p-spin = — 1 for the negative energy solutions (v, described below), and 

Np. = {Ep. + mj)^ Pi = M • ^2.15) 

Pi 

The index j denotes a quark (j = 1) or antiquark {j = 2). The values of p range from 0 to 1. 
The helicity spinors are defined in Table | for cases when the momentum is along the z axis 
(external quarks), and when the momentum is in the xz plane at an angle 6 with respect to 
the z axis. We will use a prime to distinguish the latter from the former. The v spinor, or 
negative p-spin state, used in this paper is 
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TABLE I. Helicity spinors 


external quarks 


internal quarks 




1 

2 







^ COS I 
Vsinf 



— sin I 
V cosf 


u {p,Xj) =v{-p,Xj) 


V 1 



(2.16) 


and is consistent with that nsed in Ref. |^. It is convenient to nse the helicity representation 
becanse helicity is invariant nnder rotations, and becanse the vector operator cx-p is replaced 
by scalar eigenvalnes, thns simplifying the algebra. 

The matrix element (P; is then 


0++(p, k) = (l T 4A'Apifci) , (2.17) 

where the npper sign is for the scalar vertex and the lower sign for the time-like vector case. 
We will assnme, for the time being, that the polar angle of the external qnark is 6' instead 
of 0 (as it will be later). Then 


^xy{0'0) = Sxy cos i(0 - 0') - 2A5a,-v sin 1(0-0'). (2.18) 

Therefore, forming the two independent linear combinations 

$+(p) = ^ 1^ cos ^0' - T (p, - A) sin A6*' 

<l>“(p) = T (j), a) sin A6»' -k T (^p, -A^ cos ^0' (2.19) 


Eq. (p2.11| ) becomes 

(m2+ ^ 2 ) 4>+(p) = 


(m2+ h) (p) = -0 


-0 [ ^ J!l£^v(p,k) 

J (27r)3 

|$+(/c) 1 =fTi^icos( 6'- 0') =F *f)“(/c) pifci sin(6'- 0')! 
(Ek Np^Nki 


J (27r)3 ^AEp,Ek,^^^'^^ 

|<l>“(/c) 1 =F Ti^i cos(6* — 6*') ± (k) piki sm(6 — 0')'^ . 


( 2 . 20 ) 
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Hence the interaction depends only on the difference 6 — 6' and we may set 6' = 0 without 
loss of generality. 

For later use we will record here the other p-spin matrix elements of 0 
©AA'fe k) = k) 

(2A'fci±2A:pi) =±0^+(p,fc). (2.21) 


B. Limits of the one-channel spectator eqnation 


Now we can observe that taking the limit mi —> cx) gives a Dirac equation for the light 
particle. The fixed source for the Dirac equation is a heavy quark, so the equation will 
model & Qq system, such as a D meson. As mi —>• cx), pi —0 and V (p, k) —> l/(p — k) (see 
below), giving 

(mj- /j) <I.(p) = -e / 70 i r(p - k)>l>(fc) , ( 2 . 22 ) 

where the helicity of the heavy particle was dropped because the equation is independent of 
it and we introduced the physical momentum p '2 = —p 2 = (p — Ep-^, —p) —>• (F^B,p'), with 
Eb = p — mi. In position space Eq. ( |2.22|) is 

(m 2 - + i y^i) <I)(r) = -0 D(r)<I)(r) , (2.23) 


We will return to this equation in the next subsection. 

The confining potential V (r) which appears in Dirac Eq. ( |2.23| ) is taken to be a simple 
linear potential in position space [§ 




V(r) = ar = limcrre = limcr- 

£—>0 £—>0 r 


In momentum space this potential is 


V (q) = —Sttct lim 


im < 


4e2 


0 [ (q^ + (q^ + e^Y 


(2.24) 


(2.25) 


This form of the potential is inconvenient because the limit e —> 0 must be taken numerically. 
For the Dirac equation, we use an alternative form which has the same physics 


r(q) 



(2.26) 


where q = p — k and q^ = p — k'. 
each of the terms in (^.261) is 


It is instructive to note that the position space form of 


Sna 


q" 

Sttct 


lim - 

£^0 


Sttct 


(q2 + ,2)2 

/ ^V^=lim<53(g) J 


► lim —a — 

£^o e 

Sttct 


a 

e 


ar 


+ e 


2'|2 


a 

e 


(2.27) 
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Hence the role of the delta function subtraction is to remove the inhnite constant from the 
hrst term, leaving a pure linear potential. 

In relativistic two-body equations, the potential ( |2.26|) is generalized to 0 


V(p, k) = S-ira I ^ - Ep^ 6%p - k) J 


d^k' 


Ek' (p - kr 


(2.28) 


The insertion of the energy factors is necessary to make the kernel (|2.28|) covariant, and is 
associated with the restriction of the heavy quark to its mass-shell [^. The full ICS also 
includes the covariant replacement {p — k)"^ —> {Ep^ — — (p — k)^, but in both the 

theoretical and numerical studies in this paper we have neglected retardation and use the 
simplest replacement {p — k^ —>■ — (p — k)^. We will refer to this as the quasirelativistic 
approximation, and it should be emphasized that we use this approximation in this paper 
only to simplify the discussion. We also neglect the regularization factor and form factor 
introduced in previous studies 0. 

The energy factor in the subtraction term in (|2.28| ) gives rise to a relativistic effect of 
some importance. To see this, evaluate the diagonal matrix element 


= y S^{p-k) J 

f d^r d^r 


d^k' 87rcT 


'4){k) 


(2vr) 


Ek[ {p - k'Y 


= j d\iir)[E,Y 


(2.29) 


where, in the last line, p = \/p^ a/— V^. Hence the subtraction term becomes an operator 
which is a function of p^ —> — but independent of r. It can be evaluated by standard 
means 

r d^k'Sna r d^k'Sna 

J E,, {p - k'Y = J E,, [{p - k'Y + 

a 2a i mi ml f Ep^ + p\\ 

= )} 

= - - C{pY . (2.30) 

e 

This new subtraction term contains the same singular part we found before [see Eq. ( p.27|) ] 
plus a new hnite part C{pY- The hnite part arises from the relativistic energy factor in 
Eq. ( |2.30| ), which produces an inhnitesimal modihcation of the singular part, and it vanishes 
in the mi —>■ oo limit. It has an interesting effect which will be discussed in the next 
subsection. 

The conhning potential has the property that it is very singular when q —0. This 
suggests using a peaking approximation in which 6 ~ 6', so that the coupling between 
and <I>“ can be neglected. No further approximations are needed, because we may reduce 
all factors of Ep-^ and Ek^ to derivative operators, and replace 
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(2.31) 


kiPi cos{d — 9') = 


k-p 


(Efci + mi){Ep^ + mi) 

(k-p)^ 


^2 ^2 

_ + ^ + ^ 

‘^{Eki + ^i)(-£'pi + ^i) 2 2 


because the operators k"^ (which operates on the initial wave functions) and (which op¬ 
erates on the hnal wave functions) will eventually give indentical results. Furthermore, the 
operator can be reduced to 


STTcrq^ 47ra a 

2q^ q2 r 


(2.32) 


Combining all of these effects, Eq. ( p.20| ) becomes a single Dirac-like equation with a 
momentum dependent potential 


{m2+ ^ 2 ) $(p) = 0C(p)<h(p) - 0 

1 


cPk 


(2vr: 


3 V{p, k)M^{k) 


X 


lT2(^i+P?)±2(E,,+miKE,,+mi)J ’ 


(2.33) 


where M = Np^ Nk^ / ^4Ep^ Ek^. Recalling that P 20 = M 
coordinate space form of this equation is 


Ep^ = Eb + mi - Ep ^, the 


m 2 - 


Eb + mi — Jml — 


7° -1- i'y'^di ] <F(r) 


= -0 


0 - C'(p^)] (iTpl) T 


a 


r {El + mi)" 


Ar<F(i 


(2.34) 


where we have anticipated the application to diagonal matrix elements where k'^ = p^ 


—V^, and all functions Ep^ = Ei are replaced hj \Jm\ — V^. Using Eq. (|2.34|) we can study 
the single-channel spectator equation when the mass mi is close to m 2 , and also see how it 
approaches the Dirac limit. We will study these issues approximately in Sec. III. 


C. Equations with mixed scalar and vector confinement 


Note that the operator 0 depends on its Dirac structure; it is 11 for a scalar conhnement 
and 7 *^ for vector conhnement. Hence Eq. ( p.l2|) gives 


0 = 



scalar 
vector . 


(2.35) 


In the nonrelativistic limit, the Dirac equation reduces to a Schrodinger equation for the 
upper component, and we will choose the sign of our potential so that it conhnes the positive 
energy solution in the nonrelativistic limit. Hence, in order to obtain a nonrelativistic 
conhning potential equal to err, independent of the mixing parameter y, the operator form 
(9 of a mixed kernel must be 
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(9 = (1 — y) IL ® IL — ?/7° ® 7° 


(2.36) 


Using this definition and the result Eq. (|2.34|) gives the following equation for a mixed 
confining potential 


(rn2 - [Eb + mi - Ei] 7° + iYdi^ <h(r) 

= -0 (^[ar - C(p=)] (1 - Pl(l - 2 y)) - 

Assuming a ground state solution of the form [0] 

$(r) = 


(2.37) 


fir) 

—ig{r)a ■ f 


X, 


(2.38) 


Eq. ( p.37|) reduces to the following set of coupled equations for the radial wave functions 
/(r) and g{r) 


{Eb - m 2 - [El - mi]) f + ^ + -g 

dr r 


= [{ar-C) l-pl{l- 2y) 


(T (1-2|/) 
r {El + mi)2^ 


AT/ 


{Eb + m 2 - [El - mi]) g 


dl 

dr 


= -[{ar-C) (1 - 2y) - pi 


a 


r {El + mi)2 

We will return to this coupled set of equations in the next section. 


Ug 


(2.39) 


D. Spectator equation in helicity space 


The equations we have obtained so far are convenient for approximate analysis, but 
an exact helicity decomposition is better for numerical solutions. To obtain this form of 
the one-channel spectator equation, return to Eq. (|2.9| ) and expand the wave function and 
the projection operator in terms of the helicity spinors given in Eqs. ( p.l4|) and ( p.lGj) . 
The wave function can be expanded using the decomposition of the propagator into p spin 
contributions [1^1, 


{m2+ M) 
ml - kl 


1 

2-Efc2 


E 


u{k, X2)u{k, X 2 ) 
Ek2 — ^20 — 


n(-k, A2)n(-k, A2) 
Ek2 + ^20 ~ 


Using this in Eq. (|2.8|) shows that the wave function has the form 


P^2 


(2.40) 


(2.41) 


Furthermore, the most general form of the pseudoscalar vertex function with particle 1 on 
shell is 
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u 


-(p,A)r(p)=w+(p,A) {ri7' + r27'‘(m2-^2)}, ( 2 . 42 ) 

and these Dirac operators are built only from the 2x2 matrices 11 and cx-p = 2Ap. Therefore 
in helicity space the helicity is conserved and an explicit calculation shows that 

tl+(p, A) r(p) u^ip, \,) = S,,, (2A)'« T^p) , (2.43) 

where the r^(p) are independent of the helicity. Hence the expansion (^.411) can be written 

1 


^2Ep^ T(p, A) = 


4£(p2 ^pi 


[iJia{v)u (p, A) + (2 A)^/>i6(p)m+(p, A)} , (2.44) 


where 


■ipla = 


i^lb = 


r+ 


( 2 . 45 ) 


Ek2 + Ek^ — fj, Ek2 — Ek^ + jj, 

Bringing all of these elements together, using Eq. ( |2.21|) , and choosing the sign of the 
vector interaction in accordance with Eq. (|2.36 ) gives the helicity form of the single channel 
spectator equation 

( {Eb - Ep 
\ (Eb + Ep2 


.p2 [Ep^-mi])'ilJia{p)\ ^ r y 


[Ep, - mi]) V’ib(p) ) 


Di D2\ ( 


D, 


3 D^J \ i!ib{k) j 


(2.46) 


where Eb = p — mi, 


the rescaled potential kernel is 


and Di = Ai + Bi cos 9 with 


d^k 

(^ 


Np^Np^Nk^Nk2 


and 



^P2 ^k\ ^k2 

1 

Ai = Q 

Bi = tR 

A 2 = T 2 

B 2 = ±52 

713 = 52 

B 3 = ±T 2 

A 4 = R 

B4 = tQ, 

1 +pip2hh 

R = piki 

Pj - kik2pj' 

'Ej = kj — 


(2.47) 


(2.48) 


(2.49) 


(2.50) 

PiP2kjf . 

In ( 11.49| ) the upper sign holds for scalar confinement and the lower for vector confinement 
and in (p.50|) j' ^ j. 

For the mixed scalar/vector conhnement dehned in Eq. (|2.36|) the values of Ai and 

are: 

Ai = Q Bi = -R{1 - 2y) 

A 2 = T 2 B 2 = S' 2(1 - 2y) , . 

713 = 52 B3 = T2(l-2y) ^ ^ 

A4 = R B4 = -Q(1 - 2y ). 

When the masses are equal this equation reduces to the equation previously introduced in 
Ref. 
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E. Dirac equation in helicity space 


The helicity form of the Dirac equation is obtained from ( |2.46|) by taking the mi 
limit 


oo 


f {Eb - Ep^)'ll)ia{p) \ ^ [ y (dl d 2 \( ^la{k) \ 

(Es + EpJ'0ib(p) j A Vds (i4/ \ ^/Jlbik) ) ’ 

where now V = VNp^Nj.^/ (2^JEp^Ek^ ) and dj = a* + bi cos 6 with 

ai = 1 bi = -^ 2 ^ 2(1 - 2 |/) 

02 = k 2 63 = P2(l - 2|/) 

a3=p2_ 62 = ^2(1-22/) 

a 4 , = p 2 k 2 61 = -(1-2?/). 


(2.52) 


(2.53) 


We conclude this section with a derivation of the helicity form of the Salpeter equation. 


F. Salpeter equation in helicity space 


The Salpeter m uses the approximation that the potential, or kernel, of the Bethe- 
Salpeter equation is independent of ko and po. Therefore, in coordinate space the potentials 
and the wave functions are instantaneous, i.e. ti = ^ 2 - The Salpeter Equation has two 
undesirable features. First, neglecting the energy dependence of the kernel is unphysical. 
Second, there is no Dirac limit for this equation. When the mass of one of the particles 
is taken to inhnity, the resulting equations do not reduce to a Dirac equation for the light 
quark moving in the field created by the heavy quark. Hence, it is most appropriate to use 
this equation for equal masses, far away from the one body limit. 

The direct derivation of the Salpeter equation utilizes the same steps as those used for 
the ICS equation with a few modifications. In this case pole 2, as defined in Eq. ( |2.5|) and 
Fig. P, must be included. For brevity we will only give the final result. The general equation 


IS 


rf‘«(p) = 


(Pk 

(27r)= 


K E E 0?"^ 

P 1 P 2 


(2.54) 


where pi 7 ^ p 2 and p'^ 7 ^ P 2 - The second channel wave function, denoted ip 2 a, corresponds to 
propagation of the two quarks in their negative energy state, and is equal to 

rr+ 


Ek2 + Eki + fj' 

The two wave functions, ^pla and 'ip 2 a satisfy the coupled equations 


(/? Ep2 dI!pi)Plaip) 

{p + Ep2 + Ep^)lp 2 a{p) 



Plaik) 

'<p2a{k) 


with 


Db = P1P2 + kik2 + {I - 2y) (pik2 + ^ 2 ^ 1 ) cos 6 
All other terms have the same definitions as before. 


(2.55) 


(2.56) 


(2.57) 
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III. APPROXIMATE THEORETICAL RESULTS 


In this section we develop approximations which help ns nnderstand the stability issues 
which will arise when the equations are solved numerically. 


A. Dirac solutions for large r 

We begin by studying the stability of Eq. (|2.39|) in the Dirac limit when mi —*• oo 


Eb f = {m 2 + ar) / - ^ - -9 

dr r 

EBg = {-m 2 - ar{l - 2 y)) g + 


dr 


(3,1) 


This is the exact Dirac equation for a potential which is a superposition of scalar and vector 
linear conhning forces. At large r the equations become approximately 


arf- —= 0 
dr 

-ar(l - 2 p)g+ t 0. 


(3,2) 


The solution to these equations depends on the value of y. If ?/ < 1/2, then the solution 
which approaches zero as r —> cx) is 


/(r) = Nf 5-^" g(r) = Ng , 


where 


= -Jf-2yN, 


9 ■ 


(3,3) 


(3,4) 


Note that the wave functions become less conhned as y —>■ 1/2. For y > 1/2 the solutions 
are oscillatory and escape to large r. In this case the most general solution is a linear 
combination of the following two independent solutions 


/i(r) = Nf^ sin (^ij 2 y - 1 -ar^^ gi{r) = Ng^ cos {j 2 y-l -ar^^ 

/ 2 (r) = Nf^ cos (^ij2y - 1 g 2 {r) = Ng^ sin (^ij2y-l , (3.5) 


where 


^/l- - 1 Ng^ 

Nf,= ^2y-l Ng ,. (3.6) 

This is the simple mathematical explanation behind the results shown in Figs. |] and 
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B. Estimates for the one-channel spectator equation 


Study of the solutions of the approximate one channel spectator equation, Eq. (|2.39| ), is 
complicated by the presence of the operators \Jm\ — and p\. We will therefore develop a 
variational-like method which can give us insight into the conhning behavior of the equation. 
First, if (T = 0 the exact solution of the equations is 

/W = fojoM 

9{r)=goJiM, (3.7) 


where je is the spherical Bessel function of order £, and the energy is a function of the 
parameter 7 


EbCi) = 



(3.8) 


The spectrum is continous with a gap between the positive and negative energy states. It is 
amusing to see that the energies of both the positive and negative energy states are always 
greater than the corresponding Dirac state energies, and that the negative energy spectrum 
is now bounded between —m 2 and —mi, instead of running from —m2 to — cxd. This already 
illustrates one of the new features of the ICS equation. 

When (j 7 ^ 0 we cannot solve the equation analytically, and will limit our study to the 
behavior of the expectation value of the energies as estimated by taking matrix elements 
of the equation. To compute these matrix elements we will use wave functions of the type 
shown in Fig. |^, which are constructed from spherical bessel functions of order zero and one. 
This choice makes the evaluation of functions of the operator easy. 

Ideally, the functions used should consist of a region where — is positive, and a “tail” 
region where — < 0. The functions shown in Fig. ^ were constructed from jii'yr) and 

hi{kr) joined so that the function and its first derivative are continuous. However, we found 
that the contributions from the tails did not change the qualitative behavior of the matrix 
elements, and hence we present here only the simplest results for wave functions without 
tails (where /c /7 —*• 00 , the heavy solid lines in the hgure). These results are easy to evaluate. 

Hence the “trial” wave functions we choose are 


= 


yr < TT 
yr > TT 
yr < 77-1 
yr > rii , 


(3.9) 


where je is the spherical Bessel function of order £, y is a variational parameter, and the 
constant ni = 4.493 is the location of the zero of ji. These wave functions are eigenfunctions 
of the operator V^: 


1 (9^ 

fi^)= /W 

9 {r) = - 7 ' 9 {r) • (3.10) 
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FIG. 7. The trial wave functions f{x) and g{x) as a function of x = 'jr. The different tails are 
for the cases k/'j = 0.1 (biggest tail), 0.5, and 2.0 (smallest tail), as discussed in the text. 


Hence the operators and pl can be readily calculated. 

Substituting / and g into Eq. ( p.39|) , multiplying the hrst equation by joi^r) and the 
second by jii'jr), and integrating over dPr gives the following coupled equations for /o and 
9o 


En- 


mi + 7^ — mi 


-m2^ fo + igo 


where 


and 


= (l “ Pi “ 2?/)) -^Pi(l-2r/)| M fo = Sffo 

+ 7^ -mi + m-^ go-lbfo 

= “ “ <^(7^)^ (l - ‘^y-Pl) - ^Pij -^90 = -Sggo, 


Cl = = 1.571 


C2 = 


fo x^dxj^(x) 

Jo xdxj^jx) 
Jq x‘^dxjl{x) 


c, = E = 2.669 


= 0.776 


C4 = 


x^dxjf(x) 

fo^ xdxj^(x) 
x^dxjf(x) 


= 0.412 


6 = iLfhhM = 0 734 
/;■ ■ ' 


( 3 . 11 ) 


( 3 . 12 ) 


( 3 . 13 ) 


Solving Eq. (|3.11|) gives an estimate for the eigenvalues Eb as a function of 7 , related to the 
size of the state 
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y 

FIG. 8 . The Dirac energy E = Eb as a function of the variational parameter 7 for different 
mixing ratios y = 0 (solid line), y = 0.4 (dot-dashed), y = 0.6 (dashed), and y = 1.0 (dotted). 


Eb = + 72 -mi + ^{Sf - Sg) ± ^^(2m2 + Sf + Sg^ + j'^b, (3.14) 


where Sf and Sg were defined in Eq. (|3.11|) . These energy snrfaces for a variety of cases are 
shown in Figs. S-0- In all of these cases we chose a = 0.2 GeV^ and m 2 = 0.325 GeV. We 
will now discuss some of the interesting features of these solutions. 

Note that the solutions (|3.14|) are always real, and that as 7 —0 




a 


(ci - 03(1 - 2y) ± |ci -h 03(1 - 2y)\) 


Hence the positive energy solution always approaches -|-cxo as 7 
solution goes like 


( 3 . 15 ) 

0 , but the negative energy 


E 


B 


a 


7 


-c^{l-2y) if |/< 0.795 = 


Cl + C3 
2C3 


(3.16) 


c^{2,y - 1 ) 


if 2 / > 0.795, 


and becomes positive for y > 1/2, as shown in Figs. ^ and |T^. This is a sign of instability. 
When y > 1/2 the positive energy states cannot be stable because they may always reduce 
their energy by tunneling through to a negative energy surface and sliding down to — cx). 

A similar problem may occur at large 7 , but because our estimates are less reliable here 
(we neglected the wave function tails which are more important at large 7 ) we can draw no 
hrm conclusion. As 7 —>■ cxo, 
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FIG. 9. Energies E = Eb as a function of the variational parameter 7 for the Dirac equation 
(solid line) and the ICS equation with mi = 10 m 2 (dashed line). In both cases, y = 0. 


4.0 

2.0 

E 0.0 


- 2.0 


1.0 2.0 3.0 

Y 

FIG. 10. The ICS energy E = Eb as a function of the variational parameter 7 for different 
mass ratios n = mi/m 2 = 10 (solid line), k = 5 (dot-dashed), k = 2 (dashed), and k = 1 (dotted). 
In all cases, y = 0 



21 









FIG. 11. The ICS energy E = Eb as a function of the variational parameter 7 for different 
mixing ratios y = 0 (solid line), y = 0.4 (dot-dashed), y = 0.6 (dashed), and y = 1.0 (dotted). In 
all cases ra\lm 2 = 10 . 



y 

FIG. 12. As in Fig. 11 but with mi = m 2 . 







TABLE II. Comparison of the exact and estimated solutions for the 
Dirac and ICS equations. All energies are in GeV, and the symbol - indicates 
that there is no stable solution. 


parameters 

positive energy 

negative energy 



exact 

estimate 

exact 

estimate 

mi/m 2 

V 

El 

E 

7 

E-i 

E 

7 

Dirac 

00 

0.0 

0.976 

0.950 

0.715 

-1.249 

-1.226 

0.859 


0.4 

1.028 

1.014 

0.673 

-0.660 

-0.650 

0.463 

One Channel Spectator 

10 

0.0 

0.964 

0.946 

0.635 

-1.091 

-1.034 

0.988 


0.4 

1.013 

1.007 

0.603 

-0.619 

-0.598 

0.505 

5 

0.0 

0.940 

0.926 

0.579 

-0.936 

-0.828 

1.272 


0.4 

0.992 

0.992 

0.552 

-0.548 

-0.532 

0.563 

2 

0.0 

0.857 

0.857 

0.471 

-0.607 

- 

- 


0.4 

0.928 

0.952 

0.452 

- 

- 

- 

1 

0.0 

0.745 

0.777 

0.379 

-0.330 

- 

- 


0.4 

0.853 

0.928 

0.367 

- 

- 

- 


E 


B- 


7 


(i± Vt), 


(3.17) 


and because 6 < 1 the negative energy solutions also become positive at large 7. This feature 
sets in at lower values of 7 as the mass ratio mi/m 2 decreases, as is shown in Fig. In 
fact we do note that the numerical solutions for the negative energy states are unstable for 
small values of mi/m 2 , but we see no sign of instability in the positive energy solutions for 
small values of y and all values of mi/m 2 - 

Finally, a comparison between these estimates and exact solutions for the ground state 
are summarized in Table y. Note that Eq. ( |2.39|) does a credible job of explaining the trends, 
all of which can be understood qualitatively from examination of the hgures. 

Before leaving the discussion of the ICS equation, we comment on two features of our 
estimates due to the presence of the “constant” term C{p^) of relativistic origin [recall 


Eq. (|2.30|) ]. First, note that the positive energy ICS solutions approach the Dirac limit as 
mi —>• cx) from below instead of from above, as would have been suggested by our analysis 
of the free particle case. (Note the comparison in Fig. p.) Even though the energy factor 
[El — mi] is positive, the term C is negative and is just a bit larger, giving the observed be¬ 
havior. Second, the term —C becomes more negative with decreasing mass ratio, explaining 
the drop in the binding energy as mi/m 2 decreases to unity. 
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FIG. 13. The square of the bound state mass /i^ as a function of the parameter 7 for the 
Salpeter equation with a pure scalar confining interaction (y = 0). Solid curve (m = 0.325), 
dashed curve (m = 0 . 1 ). 


C. Stability of the Salpeter equation 


Applying our technique to the Salpeter equation (|2.56|) for equal masses (mi = m2 = m) 
gives 


(/i - 2E) /i =Wifi+ W2 /2 

(ij, + 2E) f2 =-wi f2 - W2 fi (3.18) 


where 


Wi = j 

W2 = 2 


-2^p2(l-2y)| Af 

{(^ - C( 7 ^)) 2p‘(l-y) + ^p7i -2y)\M, 


(3.19) 


and we have assumed that /i and /2 are both S-states. Hence the estimated mass is 


jj? — (2E + Wi) 


Wr, 


(3.20) 


We have recovered the result that the masses always occur in ± pairs, and we see that they 
may be imaginary if \w 2 \ > \2E + tci|. 

First note that as 7 ^ 0, 


/i 


2 



(3.21) 


and as 7 —>• cx), 
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TABLE III. Comparison of the exact and esti¬ 
mated solutions for the Salpeter equation. All ener¬ 
gies are in GeV and the symbol - indicates that there 
is no stable solution. 


parameters 

m y 

exact 

El 

estimate 

7 

0.325 

0.0 

” 

0.973 

0.340 


0.4 

1.339 

1.537 

0.349 


0.6 

1.510 

1.819 

0.353 


1.0 

1.837 

2.380 

0.361 

0.650 

0.0 

3.112 

3.217 

0.466 

0.900 

0.0 

5.235 

5.396 

0.529 


/i2 = (27)2, (3.22) 

so that is always large and positive at the extreme values of 7, and must have a minimum 
for some 7. If this minimum is negative, the masses will be imaginary (i.e., the state will be 
unstable). This can occur only if m and y are small enough to satisfy the condition 

2 T1 -h I — C(7^)^ (1 — +p‘^ + 8yp^^ — (1 — 2?/)| Af 

= 2E + wi — W 2 < 0 . (3.23) 


If m = 0 this condition reduces to to 


4cr / 2\ , 

27-( Cl + C2-(1 - 2r/) < 0 . 

7 V vr/ 


(3.24) 


Hence the Salpeter equation for m = 0 is unstable only if r/ < 1/2! As m increases, this 
critical value of y decreases. If ?/ = 0, our estimate Eq. ( |3.23| ) leads to the conclusion that 
the scalar Salpeter equation is unstable only if m < 0.18; for larger values of m the equation 
has real roots for all y. This behavior is illustrated in Fig. |^, which shows that the scalar 
Salpeter equation is stable for m = 0.325 (our standard choice for the quark mass) and 
unstable for m = 0.1. 

Unfortunately, our crude estimate Eq. (|3.23|) does not reproduce the quantitative features 
of the exact Salpeter solutions as well as it did for the previous cases. The comparison 
between exact and estimated solutions is given for a few cases in Table B Note that the 
qualitative agreement is good, but that we are unable to “predict” the the critical mass at 
which the Salpeter equation becomes unstable. The exact solutions tell us that this mass is 
around 0.85 GeV, much higher than the estimated value of 0.18. 
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TABLE IV. First four positive and negative Dirac energy levels for y=0.0, 0.4, and 0.6 
with spline ranks of 20, 16, and 12. The energies are in GeV. The bold face numbers are 
unstable states with energies greater than the stable ground state, as discussed in the text. 



y = 0.0 

II 

o 

y = 0.6 

Level 

SN=20 

SN=16 

SN=12 

SN=20 

SN=16 

SN=12 

SN=20 

SN=16 

SN12 

4 

1.945 

1.945 

1.946 

2.035 

2.035 

2.035 

2.092 

2.092 

2.093 

3 

1.695 

1.695 

1.695 

1.772 

1.772 

1.772 

1.821 

1.821 

1.821 

2 

1.394 

1.393 

1.393 

1.456 

1.455 

1.455 

1.496 

1.496 

1.496 

1 

0.976 

0.976 

0.976 

1.028 

1.028 

1.028 

1.065 

1.065 

1.065 

-1 

-1.249 

-1.249 

-1.248 

-0.660 

-0.660 

-0.660 

2.028 

1.576 

1.120 

-2 

-1.575 

-1.575 

-1.574 

-0.781 

-0.781 

-0.780 

1.190 

0.861 

0.525 

-3 

-1.839 

-1.839 

-1.838 

-0.879 

-0.878 

-0.879 

0.899 

0.590 

0.278 

-4 

-2.067 

-2.067 

-2.078 

-0.963 

-0.963 

-0.964 

0.692 

0.396 

0.090 


IV. NUMERICAL RESULTS 

Now we turn our attention to the numerical solutions for the Dirac, ICS, and Salpeter 
equations. 

Numerical results are obtained by expanding the solutions in terms of splines, as de¬ 
scribed in Appendix A. In this way the integral equations in momentum space are turned 
into matrix equations and the problem reduced to a generalized matrix eigenvalue problem. 
Numerical values of the eigenvectors (expansion coefficients) and the eigenvalues (bound 
state masses or binding energies) are obtained, and the wave functions are constructed from 
the spline expansion. 


A. The Dirac equation 

The Dirac equation is reduced to the system given in Eq. ([A^ ) and Eq. ( |A^) and can be 
solved numerically on a PC in a reasonable length of time. The antiquark mass was set to 
m = 0.325 GeV and the conhnement strength a = 0.2 GeV^. We looked at four different 
values of the vector strength: y = 0.0 (pure scalar), 0.4, 0.6, and 1.0 (pure vector). The hrst 
four positive and negative energy levels for y values of 0.0, 0.4 and 0.6 are listed in Table IV 
for spline ranks of 12, 16, and 20. The pure vector case {y = 1.0) was found to be fully 
unstable, as predicted by Fig. |^, and is not listed in the Table. The eigenvalues, which for 
the Dirac equation are the binding energies, are all real and therefore pass the hrst stability 
condition (as described in Sec. ID). 

Of the four cases studied, only the negative energy levels for the y > 1/2 cases (i.e., 
y=0.Q and 1.0) vary signihcantly with the spline rank, as shown in Table IV. This violates 
the second of the stability conditions dehned in Sec. ID. Furthermore, the bold face values 
in Table IV highlight unstable states whose eigenvalues are greater than the positive ground 
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FIG. 14. Dirac positive ground state solutions for three values of the vector strength y: y = 0.0, 
£^1=0.976 GeV (circles and squares); y = 0.4, £i=1.028 GeV (solid and long dashed lines); and 
for y = 0.6, £i=1.065 GeV (heavy short dashed and dotted lines). 

state, and hence the y > 1/2 eqnations also violate the third stability condition. These 
nnstable states were identihed and tracked with changing spline nnmber by looking at their 
momentnm space strnctnre, as discnssed below. 

The Dirac wave fnnctions are shown in Figs. HM- Fig- 0 gives the positive energy 
gronnd states, Fig. 0 the hrst positive energy excited states, and Fig. 0 the negative energy 
gronnd states. By comparing the solntions for the states with y <1/2 (which are known to 
be stable) with the y = 0.6 solntions, we conclnde that (i) the positive energy y = 0.6 state 
has a strnctnre identical to the other positive energy states, and hence appears to be stable 
(as already snggested by the stability of the eigenvalne shown in Table 01), bnt (ii) the 
negative energy y = 0.6 gronnd state, shown in Fig. |^, has a radically different strnctnre 
(similar to a momentnm space delta fnnction) showing that it is indeed nnstable. All of the 
j/ = 1.0 solntions (not shown in the hgnres) have a behavior similar to the negative energy 
y = 0.6 solntion, conhrming that they are nnstable. 

The apparant stability of the y = 0.6 positive energy solntion differs from expectations 
based on the discnssion in Sec. II and examination of Fig. |^. We expect all positive energy 
solntions for j/> 1/2 to be nnstable, bnt as Fig. | shows, the positive and negative energy 
surfaces actually overlap in the y = 1.0 case but remain dearly separated for the y = 0.6 
case. This suggests that the instability of the y = 0.6 positive energy state is hard to 
observe numerically because the distance between the positive and negative energy surfaces 
is large and the “leakage” from positive to negative energy is very small (also suggested by 
Fig. §. Presumably a more precise numerical calculation would uncover some instability 
in the positive energy y = 0.6 case, but this further calculation is not needed because the 
overlap of the positive and negative energy spectrum (condition 3) is already a sign of the 
instability. 
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FIG. 15. Dirac positive first excited state solutions for y = 0.0, E 2 = 1.394 GeV (circles and 
squares), for y = 0.4, £^2=1-4:56 GeV (solid and long dashed lines), and for y = 0.6, Fi2=l-496 GeV 
(heavy short dashed and dotted lines). 



FIG. 16. Dirac negative ground state solutions for y = 0.0, E-i = —1.249 GeV (circles and 
squares), for y = 0.4, £'_i=-0.660 GeV (solid and long dashed lines), and for y = 0.6, E^i = 2.028 
GeV (heavy short dasheded and dotted lines). 








TABLE V. First four positive and negative energy levels for the ICS equation for mass 
ratios k= 5.0 and 10.0 and vector strength y=0.0 and 0.4. Here is shown in GeV and 
solutions for spline ranks of 20 and 12 are compared. 



0 

0 

II 

511 

K = 5.0 

y = 0.0 

K = 10.0 

II 

0 

K = 5.0 

0 

II 

K = 10.0 

Level 

SN=20 

SN=12 

SN=20 

SN=12 

SN=20 

SN=12 

SN=20 

SN12 

4 

2.109 

2.113 

2.073 

2.078 

2.225 

2.227 

2.165 

2.168 

3 

1.808 

1.808 

1.783 

1.783 

1.898 

1.899 

1.858 

1.858 

2 

1.443 

1.443 

1.435 

1.435 

1.509 

1.509 

1.495 

1.494 

1 

0.940 

0.939 

0.964 

0.964 

0.992 

0.992 

1.013 

1.013 

-1 

-0.936 

-0.936 

-1.091 

-1.090 

-0.548 

-0.569 

-0.619 

-0.619 

-2 

-1.084 

-1.084 

-1.333 

-1.332 

-0.570 

-0.607 

-0.715 

-0.715 

-3 

-1.173 

-1.170 

-1.511 

-1.515 

-0.600 

-0.637 

-0.786 

-0.785 

-4 

-1.233 

-1.259 

-1.650 

-1.642 

-0.630 

-0.675 

-0.841 

-0.848 


We conclude that the fourth stability condition largely reinforces the conclusions we 
have already drawn, but that it should be used in conjunction with the other three. The 
stability of a single state cannot easily be determined solely by tracking (with changing 
spline number) its behavior. A reliable conclusion requires the examination of the entire 
spectrum, with particular attention to condition 3. 


B. The one-channel spectator eqnation 

As in the Dirac case the antiquark mass will be set to m 2 = 0.325 GeV and the conhne- 
ment strength to a = 0.2 GeV^. We will present results for heavy quark masses mi = nm 2 
with the mass ratio k, = 10, 5, 2, and 1. In order to compare the IGS results to those 
obtained from the Dirac equation, we dehne an effective Dirac-like binding energy Eo using 
the relation 


fi = Ed + mi 


(4.1) 


where fi is the mass eigenvalue obtained from the IGS equation. This relation insures that 
the effective IGS binding energy must approach the Dirac binding energy as mi —> 00 . 
Tables V and VI give these effective binding energies (instead of the bound state masses). 

Note that results for the equal mass (k = 1.0) IGS equation are included only for 
comparison because the IGS equation should not be used for equal mass systems. If the 
equal mass particles are identical (as in NN scattering) the equation must be symmetrized 
in order to preserve the Pauli principle. Even if the equal mass particles are not identical, 
as for the qq pairs discussed in this paper, the equation must still be symmetrized to insure 
charge conjugation invariance. Furthermore, for bound states with a very small mass (e.g., 
the pion) the symmetrized two channel spectator equation defined in Ref. 
used. 


should be 
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TABLE VI. First four positive and negative energy levels for the ICS equation for the 
mass ratio k= 1.0 and vector strength y=0.0 and 0.4. Here is shown in GeV and solutions 
for spline ranks of 24, 20, 16, and 12 are compared. 




y = 

0.0 



y = 

0.4 


Level 

SN=24 

SN=20 

SN=16 

SN=12 

SN=24 

SN=20 

SN=16 

SN12 

4 

1.881 

1.881 

1.881 

1.881 

2.222 

2.222 

2.222 

2.223 

3 

1.630 

1.630 

1.630 

1.632 

1.884 

1.884 

1.884 

1.884 

2 

1.294 

1.293 

1.293 

1.293 

1.461 

1.461 

1.461 

1.461 

1 

0.745 

0.745 

0.745 

0.745 

0.853 

0.853 

0.853 

0.853 

-1 

-0.329 

-0.330 

-0.331 

-0.334 

0.933 

0.724 

0.508 

0.284 

-2 

-0.331 

-0.332 

-0.335 

-0.341 

0.727 

0.527 

0.326 

0.122 

-3 

-0.334 

-0.337 

-0.342 

-0.354 

0.577 

0.387 

0.196 

0.005 

-4 

-0.338 

-0.343 

-0.353 

-0.379 

0.454 

0.272 

0.091 

-0.087 


The eigenvalues are real for all values of the vector strength y and the mass ratio k 
(condition 1). However only systems with a vector strength less than 1/2 (0.0 and 0.4) have 
stable eigenvalues (condition 2). Cases which fail the hrst two stability conditions {y = 0.6 
and 1.0) are not listed in the eigenvalue tables. Table 0 shows the eigenvalues for mass 
ratios fi:=5.0 and 10.0. These cases are very similar to the Dirac cases, and the table shows 
that in all cases the spectra satisfy condition 3 (no overlap of the positive and negative 
energy sectors). Table |VT| shows the eigenvalues for the equal mass case (k= 1.0). Note that 
condition 3 is violated for |/=0.4; at a spline rank of 24 the negative energy state (shown in 
bold face) crosses into the positive energy sector. In the equal mass case only the pure scalar 
interaction is stable. The binding energies for k = 2.0 (not shown in the tables) exhibit the 
same behavior as for k = 1.0. 

Wave functions for the ICS equation are shown in Figs. 17-^3]. In Figs. P!7HTU| the wave 
functions for a large mass ratio and a pure scalar conhnement are compared with the Dirac 
solutions. Both the positive and negative states for these systems are completely stable and 
very similar to the corresponding Dirac solutions. We also observe how the ICS binding 
energies approach the Dirac values as k is increased. 

Figs. 0 and 0 show the wave functions for large mass ratios and a vector strength of 
0.4. For fi:=10.0 the system is once again totally stable, while for fi:=5.0 only the positive 
states are stable. In this case the instability of the negative energy states is not accompanied 
by violation of condition 3; the only indication of instability is the variation of the negative 
energy levels with spline rank (condition 2), as shown in Table 0. In this case the structure 
(condition 4) reinforces condition 2, and we have a hrst example of a system where the 
positive energy solutions are stable and the negative energy ones are not. 

The positive and negative ground states for k = 1.0 and 2.0 are shown in Figs. and 
Note that the positive energy states are stable while the negative energy ones are not. 
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Here the instability of the negative energy states is only apparant from an examination of 
the structure of the wave functions; neither condition 2 (variation of the energy with spline 
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p (GeV) 

FIG. 17. Positive ground state soutions for the quasirelativistic ICS equation with a pure 
scalar interaction. The solid and long dashed lines are for k=5.0, 771 =0.940 GeV; the heavy short 
dashed and dotted lines are for k= 10.0, i7i=0.964 GeV. The scalar ground state Dirac solution for 
El = 0.976 GeV is shown for comparison (circles and squares). 



p (GeV) 

FIG. 18. Positive first excited state solutions labeled as in Fig 17. Here the k = 5.0 solution has 
an energy of E 2 = 1.443 GeV and the k = 10.0 solution an energy of E 2 = 1.435 GeV compared 
to the Dirac energy of E 2 = 1.394 GeV. 
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FIG. 19. Negative ground state soutions labeled as in Fig 17. Here the k = 5.0 solution has an 
energy of E^i = —0.936 GeV and the k, = 10.0 solution an energy of E^i = —1.091 GeV compared 
to the Dirac energy of E^i = —1.249 GeV. 



FIG. 20. Positive ground state solutions of the quasirelativistic ICS equation with a mixed 
scalar and vector interaction {y = 0.4) for two mass ratios k. The solid and long dashed lines 
are for k = 5.0, Ei = 0.992 GeV, and the heavy short dashed and dotted lines are for k = 10.0, 
El = 1.013 GeV. The circles and squares show the solution for the Dirac equation with Ei = 1.028 
GeV. 
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FIG. 21. Negative ground state solutions of the quasirelativistic ICS equation for y = 0.4 
labeled as in previous figure. Here k = 5.0, E^i = —0.548 GeV and k = 10.0, E^i = —0.619 GeV. 
The comparison Dirac level has energy E-i = 0.660 GeV. 



p (GeV) 

FIG. 22. Positive ground state soutions for the quasirelativistic ICS equation with a pure 
scalar interaction. The solid and long dashed lines are for k=1.0, £'1 =0.745 GeV; the heavy short 
dashed and dotted lines are for /t=2.0, £1=0.857 GeV. The scalar ground state Dirac solution for 
£i = 0.976 GeV is shown for comparison (circles and squares). 
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p (GeV) 

FIG. 23. Negative ground state solutions labeled as in previous figure. Here the k = 1.0 solution 
has an energy of E-i = —0.330 GeV and the k, = 2.0 solution an energy of E-i = —0.607 GeV 
compared to the Dirac energy of E-i = —1.249 GeV. 

rank) nor condition 3 (penetration of the positive energy sector) seems to occnr. 

In conclusion, the ICS system becomes more stable as the vector strength is decreased 
and the mass of the heavy quark is increased. This will be summarized further at the end 
of this section. 


C. The Salpeter Equation 


The use of pure scalar conhnement with the Salpeter equation gives the hrst example 
of instability due to the mass eigenvalues becoming complex (condition 1). Actually, the 
eigenvalues become pure imaginary because the mass squared is real and negative. This 
situation is accompanied by a very rapid variation of with spline rank, as shown in Table 
Vlf| . However, for y=0A the tabulated spectra do not vary with the spline rank, and these 
states are stable, as shown in Figs. |2^-^. Fig. ^ also shows that the wave functions 
for positive and negative energies are identical provided 'ijjia ^ '02a- This is a further 
consequence of the symmetry of the Salpeter equation which produces pairs of eigenvalues 
with the same magnitude and opposite signs. 

The two hgures. Fig. ^ (ground state) and Fig. |2^ (second excited state), demonstrate 
that these Salpeter systems have solutions comparable to their Dirac counterparts. In addi¬ 
tion, Fig. ^ illustrates that the y=0.6 and 1.0 solutions are indeed stable by showing that 
they have the correct structure with the right number of nodes for a second excited state. 

While it is true that the scalar Salpeter equation is unstable for equal quark-antiquark 
masses of 0.325 GeV, increasing the mass of the quarks will give stable solutions (this was 
anticipated by the discussion in Sec. III). We hnd that the lower mass states of the |/=0.0 
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TABLE VII. Square of the mass in GeV^) of the first four levels of the Salpeter 
equation for y=0.0 and 0.4 with various spline ranks. 


Level 

o 

o 

II 

II 

o 

SN=20 

SN=12 

SN=20 

SN=16 

SN=12 

4 

0.685 

2.173 

5.632 

5.632 

5.674 

3 

-1.074 

1.538 

4.383 

4.383 

4.385 

2 

-3.869 

0.931 

2.977 

2.977 

2.976 

1 

-8.705 

-0.051 

1.339 

1.339 

1.339 



p (GeV) 


FIG. 24. Positive and negative ground state solutions for the y = 0.4 quasirelativistic equal 
mass Salpeter equation, /xi=1.157 GeV (solid and long dashed lines) and /i_i=-1.157 GeV (heavy 
short dashed and dotted lines). The positive ground state Dirac solutions for y=0A, i7i=1.028 
GeV (circles and squares) are shown for comparison. 
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FIG. 25. The positive second excited state solutions for the y=0.4 equal mass Salpeter equation, 
/is = 2.094 GeV (solid and long dashed lines) are compared to the second positive excited state 
Dirac solution for y = 0.4, £'3=1.772 GeV (circles and squares). 



p (GeV) 

FIG. 26. Positive second excited state solutions for the Salpeter equation for a variety of 
scalar/vector mixings: pure vector y=l.O, /i3=2.565 GeV (circles and squares); y=0.G, /i3=2.284 
GeV (solid and long dashed lines); and y=0A /i3=2.094 GeV (heavy short dashed and dotted lines). 


36 













FIG. 27. Stable ground state solution for the Salpeter equation with a pure scalar confining 
interaction. In these cases m = 0.85 GeV and the Salpeter bound state has mass /ii=2.185 
GeV. The solutions for SN=20 (circles and squares) and SN=30 (solid and long dashed lines) are 
compared. (Note that the solutions are not normalized - see the discussion in the text.) 

Salpeter equation are stable when the quark mass is increased to m=0.85. The ground state 
wave functions for this case are shown in Fig. 0. where solutions for spline ranks of 20 and 
30 are compared (since the wave functions have not been normalized, only the shape of the 
two solutions should be compared). Solutions obtained for somewhat lower masses (m=0.65, 
for example) appear stable for SN=20, but the spectrum shows some instability for SN=30. 
In general, the number of stable states for the pure scalar Salpeter equation increases as the 
quark mass increases. Further study is needed to obtain a detailed understanding of the 
stability of the purely scalar Salpeter equation. 


V. CONCLUSIONS 


Table |V111| summarizes the results presented in the previous sections, which are also 
outlined below: 


• The Dirac equation is stable if the scalar conhnement is stronger than the vector 
conhnement {y < 1/2). 

• The Salpeter equation is stable if the interaction is mostly vector, and perhaps also for 
pure scalar exchanges with a large quark mass. The precise boundary between stable 
and unstable solutions is presumably a function of the quark mass m and the vector 
strength y, and we have not mapped it out. 

• The one channel spectator (ICS) equation has the Dirac limit, as expected. This means 
that for large mass ratios k = mi/m 2 , it is stable if the interaction is predominately 
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TABLE VIII. Stability Results (the table lists the region of stability or the first of the 
four tests that the system fails). 



o 

o 

II 

y = 0.4 

y = 0.6 

2/= 1.0 

Dirac 

stable 

stable 

Cond 2 

Cond 2 

ICS K=1.0 

positive 

Cond 3 

Cond 2 

Cond 2 

ICS k=2.0 

positive 

Cond 3 

Cond 2 

Cond 2 

ICS k=5.0 

stable 

positive 

Cond 2 

Cond 2 

ICS K=10.0 

stable 

stable 

Cond 2 

Cond 2 

Salpeter 

m> 0.85 GeV 

stable 

stable 

stable 


scalar {y < 1/2). However, as the mass ratio decreases toward unity, the region of 
instability grows. As we decrease k for a fixed vector strength y < 1/2, the negative 
energy states will hrst become unstable, and then the positive energy states may follow. 
However, if the vector strength is small enough (e.g., y = 0) the positive energy states 
appear to be stable for all mass ratios. 

The usefulness of an equation where only part if the spectrum is stable depends on whether 
or not the spectrum of unstable states is clearly separate from the spectrum of stable states 
(i.e. Condition 3 is met). The ICS equation for scalar conhnement has this feature; the 
unstable states are those which map, in the Dirac limit, into negative energy states. If one 
is content to exclude these states from consideration on physical grounds then the scalar 
ICS equation can be used to describe conhned Qq systems for all mass ratios. The Salpeter 
equation can also be used for equal mass qq systems unless the conhnement is predominately 
scalar and the quark masses are not large. 

This conclusion answers one of the questions raised in the introduction; clearly the 
stability of vector or scalar conhnement depends on the relativistic equation used. Scalar 
conhnement is stable if the ICS equation is used and vector conhnement is stable if the 
Salpeter equation is used. 

We emphasize that our study of the stability of the spectator equation is preliminary for 
three reasons: 

• Only the ICS equation has been studied. As emphasized before, a two channel spec¬ 
tator equation must be used if the bound state mass is small (the pion), and any 
spectator equation must be explicitly symmetrized if the quark masses are equal. 

• Our study of the ICS equation was limited to the quasirelativistic approximation, in 
which retardation is neglected. However, neglecting retardation usually leads directly 
to the Salpeter equation, and the attempt to include it (at least approximately) is the 
principle reason for choosing to use a spectator equation in the hrst place. Including 
retardation in our analysis (planned for a later work) may alter our conclusions. 
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• Only the time-like part of a vector confinement (i.e. has been studied. There are 

preliminary indications that our results will change when the full vector interaction 
7 ^ 7 ^ is included. 


Our results for the Salpeter equation agree with Ref. Q], but disagree with the results 
obtained by Parramore and Piekarewicz 0, who found the Salpeter equation to be unstable 
once the vector strength dropped below one-half, regardless of quark mass. However, as 
stated above, we hnd that the Salpeter equation is stable for a vector strength 0.4, and is 
even stable for a pure scalar interaction provided the quark mass is sufficiently large. We 
looked at one of the cases they found to be unstable (<7=0.29, m=0.9 GeV, with 25 basis 
states), and found it to contain stable states. A possible explanation for this difference is 
that we use cubic splines for our basis functions, while nonrelativistic harmonic oscillator 
wave functions were used in Ref. 0. 

There are other equations which can be used to model the quark-antiquark system. 
Tiemeijer and Tjon |]^ explored two such equations, the Blankenbecler-Sugar-Logunov- 
Tavkhelidze (BSLT) ||T^ equation and the equal-time (ET) equation of Wallace and Man- 
delzweig |^. The kernels for both equations contained one-gluon-exchange (with the full 


four vector structure) and a linear conhning term (with a mixed scalar-four vector struc¬ 
ture). They found that increasing the vector strength of the conhning term improved the 
phenomenology, but that some mesons became unstable for vector strengths of more than 
about 0.25, depending on the equation and gauge used. These results reinforce the general 
conclusions of this paper: stability depends on both the Lorentz structure of conhnement 
and on the type of relativistic equation used. 

We have seen that the study of the mathematical stability of relativistic equations re¬ 
quires the examination of both local and global features of the eigenvalue spectrum and have 
introduced four conditions which must be satished for an equation to give stable solutions. 
Using these stability criteria we hnd that the Lorentz structure of the kernel and the equa¬ 
tion used to model the meson both play a crucial role in the mathematical stability of the 
system. Clearly further research on this topic is needed. 
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APPENDIX A: SPLINE FUNCTIONS AND NUMERICAL METHODS 


To solve the equations in this paper numerically, we expand each momentum space wave 
function in terms of cubic splines 

SN 

V'i(p) = > (Al) 
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where a* are the expansion coefficients (which become the eignevectors of the problem), jSj 
are the spline fnnctions, and SN is the nnmber of spline fnnctions in the expansion (the 
spline rank). In all of the equations studied there are only two independent wave functions, 
so i = 1 or 2. Since the angular integrations are performed analytically, the wave functions 
depend only on the magnitude of the momentum p. Once Eq. ([Al| ) is substituted for each 
of the wave functions, both sides of the equation are operated on by the integral operator 

jp^f3i{p)dp. (A2) 


This reduces the integral equations to matrix equations with dimension 2SNx2SN and of 
the general form 


A, 


ij 


0 


0 A 


ij 


+ 


Bl! 


R22 

Blj 






y^f yf 



(A3) 


These equations are then solved for the eigenvalues A and the eigenvectors {aj, a|}. In the 
following subsections we give the forms of the matrices A and V for each case studied in 
this paper. 


1. Dirac equation 


The Dirac equation was given 
functions are 


in Eq. (|2.52|) and (|2.26|) . 


The two independent wave 


V’l = V'la 
^2 = -016 , 


and A = Eb, 

roo 

Aij= P^dp j3i{p)pj{p), 

Jo 

and 

poo 

B]} = -B’ff = - P^dp Bp(di{p)Pj{p ). 

Setting m 2 = m and using the notation 


N 

Mp) = -?= (di{p), 

2,Ep 

the potential matrix can be written 



L L ixl Cs) ■ 


(A4) 

(AS) 

(A6) 


(A7) 


(AS) 
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The functions rj and C, are 


r]i = ai + bi Ci = bi 


(A9) 


where a* and bi were dehned in Eq. and if rji = 7]i{p,k), then rj'- = rii{p,p). The 

fnnctions Vq and Vi are 


k) = - f" dz 

’ 2 7-1 (p2 _ 2pkzY {p^ + k‘^Y ~ ‘^p‘^k‘^ 

1 /-i [z-l) _ 1 pk 

VPP, K) 2 ^p2 + ^2 _ 2pkzy 2 (j 92 + F + 2j9fc) 

1 . (p^ + k‘^ + 2pk \ 
3 ^^2 _|_ /j2 _ J 


(AlO) 


2. One-channel spectator equation 


The ICS equation in helicity form was given in Eq. (|2.46|) with the potential dehned in 
Eq. (|2.28|) [with (p — ky —^ (p — k)^ as discussed in Sec. IIB]. The two independent wave 
functions are as in Eq. (A4|) and A = /i = mi + Eb- The matrix A is identical to the Dirac 
case, but now 


POO 

^ij = ~ / A(p)/?i(p) 

*7 0 

POO 

B‘if= - / p^dp (Epj - EpJ I3i{p)l3j{p ). 

*7 0 

Introducing the notation 


Fi{p) 


Api Ap2 

Ep^ 


A(p), 


(All) 


(A12) 


the potential matrix can be written 


with 


K 


4cr / 

TT 

4cr r 
TT 7o 


0 




dpdkV„{p,k)F,{p)I^F,{k)(^'l^^ 

V'J 

1 

dpdkVi{p,k)Ei{p)Ej{k)Q^^ , 


(A13) 

Vi = Ai + Bi Ci = Bi, 


(A14) 


where the Aj and i?j were dehned in Eq. (|2.51|) . The meaning of the prime in p' is the same 
as in p' and Vq and Vi are as before. 
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3. Salpeter equation 


The Salpeter equation is given in Eq. (|2.56|), with the masses both equal to m. Now 


= -Ipla 
'lp2 = 'lp2a , 


(A15) 


and X = fi. The matrix A is identical, but B is two times the Dirac B. The potential matrix 
is similar to Eq. (|A13|) with 


V2 ( 2 ^ Cs 

Vs -V5 Cs ^ -C 5 

Vi -Vi Ci -Cl, 

and, from Eq. ( p.57|) , 

175 = + P + 2(1 - 2y)pk Cs = 2(1 - 2y)pk . 


(A16) 


(A17) 


4. Splines 


The solution to the wavefunctions used in this paper are based on a set of third order 
polynomial functions called cubic splines. Used previously in papers such as Ref. [|^], they 
have proven versatile enough to model all of the wavefunctions examined in this paper. 

The wave function expansion was given in Eq. (|A1D - Each spline is constucted from four 
separate functions. The function used depends on the argument and the spline index j as 
shown: 


4(3j{x) 


' {x-Xj-2p 
^ ’ 

^ I Q(x-a;,_i) I ^(x-Xj-iP 

^ h , 

< 1 I Q (Xj + I-X) Q {X^ + I-Xy _ Q {X^ + I-Xf 

{Xj + 2-xf 
1? ’ 

0 


X e [Xj-2,Xj-i] 
X G [xj-i,Xj] 

X G [xj,Xj+i\ 

X G [Xj+i,Xj+2] 
otherwise. 


(AI 8 ) 


Each spline is defined on the interval from zero to one. This range is divided into sectors 
whose size, h = 1/(SN +1), depends on the spline rank. Each sector is bounded by nodes at 
Xk and Xk+i, with the number of nodes equal to SN+2. The hrst node, Xi, is always located 
at zero, and the last one, xsn +25 at one. The spline curves for a spline rank of 4 are given in 
Fig. 1^. The standard choice for our calculation was a rank of 20 (20 splines in each wave 
function expansion). 

None of the nodes may he outside of the interval from 0 to 1, so the first spline, j = 1, 
is defined entirely by the third and forth functions given in Eq. (|A18|) . It has a zero slope at 
X = 0. The j = 2 spline was defined in a special way so that it too will have zero slope at 
X = 0 (insuring that all the splines have this property). To accomplish this the first sector 
(which lies between xq and Xi and is hence outside the acceptable range of support) will be 
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FIG. 28. Spline rank 4 curves (1 solid, 2 long dashed, 3 short dashed, and 4 dotted). 


“folded over” onto the interval between [xi,x- 2 \. Hence, in the interval between [xi,X 2 ] the 
second spline is defined to be 


4A(x) = 1+3^:;^+ 


h 








(A19) 


This is an exceptional case, and all other splines are dehned following Eq. ( A18 ) in a straight¬ 
forward fashion. 

The splines defined in Eq. ( A18|) are only continuous up to their second derivative. 
Therefore, in order to obtain convergence the integrals must be separately evaluated for 
each sector, and the results from all the sectors summed up afterwards. Special care must 
be taken in evaluating those contributions to the double integral of the potential which 
include singularities. These are evaluated by choosing points equally spaced on each side of 
the singularity so that a well defined limit is obtained. 

To use the splines to describe the wave functions, the interval [0, cx)) is mapped into the 
line segment [0,1] using the tangent mapping 


x=2arctanf|: 
TT VA 


(A20) 


with A = 1 GeV. This mapping alters the shape of the splines, as illustrated in Fig. 29 


When the spline rank is increased the sectors become smaller and the range in momentum 
space over which the splines are signihcantly different from zero increases. Thus, the wave 
function is more accurately modeled as the spline rank increases. Of course this higher 
precision must be balanced by consideration of computation time. 
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FIG. 29. Spline rank 4 curves (1 solid, 2 long dashed, 3 short dashed, and 4 dotted) with 
momentum arguement. 
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